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Abstract. The two-dimensional n-body problem of classical mechanics is a non- 
integrable Hamiltonian system for n > 2. Traditional numerical integration algorithms, 
which are polynomials in the time step, typically lead to systematic drifts in the 
computed value of the total energy and angular momentum. Even symplectic 
integration schemes exactly conserve only an approximate Hamiltonian. We present 
an algorithm that conserves the true Hamiltonian and the total angular momentum to 
machine precision. It is derived by applying conventional discretizations in a new space 
obtained by transformation of the dependent variables. We develop the method first 
for the restricted circular three-body problem, then for the general two-dimensional 
three-body problem, and finally for the planar n-body problem. Jacobi coordinates are 
used to reduce the two-dimensional n-body problem to an (n — l)-body problem that 
incorporates the constant linear momentum and center of mass constraints. For a four- 
body choreography, we find that a larger time step can be used with our conservative 
algorithm than with symplectic and conventional integrators. 
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1. Introduction 

The n-body problem is the study of the motion of n arbitrary particles in space according 
to the Newtonian law of gravitation. When n = 2 (the Kepler problem), the problem 
has a well-known analytic solution, but Poincare has shown that the system is in general 
non-integrable for n > 2. To approximately solve these cases, one often attempts to 
discretize the equations of motion and study the evolution of the system numerically. 
However, discretization of a system of differential equations typically leads to a loss of 
accuracy; first integrals of the motion may no longer be preserved and the phase portrait 
may become inaccurate. This often necessitates the use of small time steps, so that many 
iterations will be required. In this article, we demonstrate that conservative integration 
can be used to obtain an accurate picture of the dynamics even with a relatively large 
time step. 

Conservative integration was introduced by Shadwick, Bowman, and Morrison 
IjH, 0, |H . These authors argued that a more robust and faithful evolution of the dynamics 
can be obtained by explicitly building in knowledge of the analytical structure of the 
equations; in this case, by preserving the known first integrals of the motion. They 
illustrated the method applied to a three-wave truncation of the Euler equations, the 
Lotka-Volterra problem, and the Kepler problem. In this work, we extend the method 
to the equations of motion of n bodies in space, first to the circular restricted three-body 
problem, then to the general three-body problem, and finally to the full n-body case. 
For simplicity we only consider two-dimensional motion (a reasonable assumption for 
all of the planets in the solar system except for Pluto); extending this work to three 
dimensions should be straightforward. 

2. Conservative Integration 

The equations describing the motion of the solar system form a conservative system: the 
friction that heavenly bodies sustain is so small that virtually no energy is lost. Both 
the total energy and total angular momentum are conserved. We argue that a robust 
integration algorithm should preserve both of these invariants. 

One way to accomplish this is to transform the dependent variables to a new space 
where the energy and other conserved quantities are linear functions of the transformed 
variables, apply a traditional integration algorithm in this space, and then transform 
back to get new values for each variable |2j. This approach is motivated by the 
following trivial lemma. 

Lemma 1 Let x and c be vectors in~R n . If f : M n+1 — > M n has values orthogonal to c, so 
that I = cx is a linear invariant of the first- order differential equation dx/dt = f(x,t), 
then each stage of the explicit m-stage discretization 




j = 1, . . . , m, 



(1) 



k=0 
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also conserves I, where r is the time step and bjk € M. 

Proof. For j — 1, . . . , m, we have c • Xj = c-xq + t Yjk=o ^jk c ■ f( x k, t + a^r) = c • a?o- ° 
A conservative integration algorithm can be constructed by writing any 
conventional integration algorithm of the form ([!]), for which specific values of aj and bjk 
are known, in a transformed space. For example, consider the second-order predictor- 
corrector (2-stage) scheme for evolving the system of ordinary differential equations 
dx/dt = f(x, t), 

x = x + rf(x ,t), (2a) 

x(t + r)=x + £[/(aj , t) + f(x, t + r)], (2b) 

where we now write x instead of X\. In the conservative predictor- corrector algorithm, 
one seeks a transformation £ = T(x) of the dependent variable x such that the quantities 
to be conserved can be expressed as linear functions of the new variables £j, i = 1, . . . , n. 
Then, keeping (pc|) as the predictor, in the transformed space one applies the corrector 

£(t + t) = £ + ^[T'(x )f(x Q , t) + T'(x)f(x, t + r)}, (3) 

where £o — T(x ) and T' is the derivative of T. The new value of x is obtained by 
inverse transformation, x(t + r) = T _1 (^(t + r)). Often the inverse transformation 
involves radicals, and if the argument of the radical becomes negative, it is possible to 
use a finite number of time-step reductions to integrate the system through this region 
@; this approach is particularly advantageous when the time step is chosen adaptively. 
Another way to deal with noninvertible transformations is to switch to a conventional 
(e.g. predictor-corrector) integrator for that one time step. If the inverse transformation 
involves several branches (e.g. because of a square root), the correct branch can be 
distinguished with sufficient accuracy using the conventional predictor solution. The 
error analysis for a second-order predictor-corrector algorithm is described in [Appendix 
[A[ Higher-order conservative integration algorithms are readily obtained in the same 
way, by coupling the first m — 1 "predictor" stages from ([!]) with the final conservative 
corrector stage 

m—1 

£{t + t) = £ + r b mk T'(x k )f(x k , t + a jT ). (4) 

k=0 

According to Iserles @j, a major drawback of traditional non-conservative 
integration is that numbers are often "thrown into the computer." Mathematical models 
are often discretized according to algorithms that have little to do with the original 
problem. Iserles argued that one should develop computational algorithms that reflect 
known structural features of the problem under consideration (e.g. see 0, ||). The 
conservative predictor-corrector is an example of such an integrator. In the examples 
given by [1], |3], the transformation T is tailored to the system at hand; there is obviously 
no generic transformation that can be used to integrate an arbitrary conservative system. 
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It is interesting to compare conservative integration (which conserves the value of 
the Hamiltonian) with symplectic integration (which conserves phase-space volume; see 
Refs. ||7|1, p, flOf, and According to Ge and Marsden (1988), if an integrator 



is both symplectic and conservative, it must be exact. Normally we do not have the 
luxury of an exact discretization at our disposal. The drawback then with conservative 
integration is that the Hamiltonian phase-space structure will not be preserved, just 
as for symplectic integration the total energy will not be conserved. Which method is 
preferable depends on the physical structure of the problem being investigated. 

Another important advantage of conservative integration algorithms is that, unlike 
typical symplectic integration schemes, they are explicit. Although in some cases the 
inverse of the transformation T may be defined by an implicit equation that requires 
iteration to solve (using the predicted value as an accurate initial guess), this is really 
nothing more than a special function evaluation; the time-stepping scheme itself, being 
causal, is explicit. 

With conservative integration, one can preserve all of the known invariants of the 
n-body problem conserved exactly, even for large time steps. This can lead to a more 
accurate picture of the motion of the bodies |], figure 9] for the same computational 
effort. In the next section, we motivate the extension of the method of conservative 
integration to the n-body problem by briefly revisiting the treatment of the Kepler 
problem in Ref. 

3. Kepler Problem 

The Kepler problem describes the motion of two bodies with masses mi and m 2 located 
at positions r\ and r 2 , respectively. The dynamics can be reduced to an equivalent one- 
body problem, the behaviour of a single particle of mass m = mim 2 / '(mi + m 2 ) at the 
position r = r 2 — r\ under the influence of a central gravitational force. This force may 
be expressed as the gradient of the potential function V = —k/r, where k = Gm 1 m 2 
and G is the universal gravitational constant. The equations of motion can be written 
in terms of the radial velocity v r and the polar coordinate angle 9 of the particle, 

— = v r , (5a) 
dv r f 1 fdV\ 

(56) 

„ ■>■ ( 5c ) 

at mr z 

where I is the (constant) total angular momentum. It is convenient to rewrite the 
equations in terms of the linear momentum p = mv r and the angular momentum t. 

dr dH p N 



dt m 2 r 3 m \ dr 
d6 I 
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dt d£ mr 2 ' 

di dH . JX 



where the Hamiltonian 

is also conserved. 
3.1. Integration 

To set the framework for extending two-body conservative integrators to the n-body 
problem, we slightly generalize the presentation in Ref. |3J to make the constant I a 
variable that is formally integrated, but which remains constant. 

The predictor step of the conservative integrator is given by (|2o|), where x = 
(r, 9,p, £). To derive the corrector, the vector (r,p, £) is transformed to (£i, £2, £3), where 

6 = --, (8a) 
r 

p 2 i 2 

^ 2 2m 2mr 2 ' ^ ^ 

6 = i- (8c) 

On differentiating these equations with respect to time and exploiting the fact that both 
H = £1 + £2 and L = £ 3 are conserved, one finds 

6 = % (9a) 
mr A 

6 = -L (96) 
£3 = 0. (9c) 
After applying ([3]), the inverse transformation 

r = ~, (10a) 

£ = 6, (106) 



£ 2 

p = sgn(p) d 2m£ 2 - (10c) 

is used to update the values of the original variables at the new step. See Ref. |]J for 
details on how the invariance of the Runge-Lenz vector A = v x £ + Vr is exploited to 
evolve #.| 

Before generalizing the integrator of Shadwick et al. to the n-body problem, it is 
instructive to consider first the special case of the restricted three-body problem. 
J There is a typographical error in Eq. (54b) of Ref. ; it should read 



v r (t + r) = sgnK) Jt* + - 2 [y 2 - -^—^ j - 2 -. (11) 
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4. Restricted Three-Body Problem 

Suppose that two bodies of masses m\ and mi, called the primaries, revolve around their 
center of mass in circular orbits. The circular restricted three-body problem describes the 
motion of a third body, with a mass 772.3 that is negligible compared to mi and 777.2, at 
coordinates (x, y) in the plane of motion of the other two bodies. The third body does 
not influence the motion of the other two. The derivation of the equations of motion 
for the restricted problem is described in fl3| . The Hamiltonian is given by 

H = \ (A 2 + y 2 ) - \ (y 2 + x 2 ) - ^ - A (12) 

where r\ = (x — fi) 2 + y 2 , r\ = (x + 1 — fi) 2 + y 2 , and ji = rrtij (mi + m 2 ). In terms of 
the canonical variables 

qi = x, q2 = y, Pi = x-y, p 2 = y + x, (13) 
the Hamiltonian appears as 

H = n(Pl +^2) +PlQ2 ~P2?1 - --— . (14) 

2 r 1 r 2 

The equations of motion are then 

Qi = tt- = Pi + 92, (15a) 

OH . . 

92 = ^— =P2-<h, (156) 

Op 2 



dH 1 — fi f N 

<9gi 



Pi = ~— =^2 3— (?1 - -3(91 + 1 ( 15c ) 



dH li 

P2 = -TT- = ~Pl 3—92 - — <?2, (15d) 

0*^2 ?i 
and the Hamiltonian can be rewritten as 

^.i. Integration 

The conventional predictor for this system is 

Qi = Qi + ftT, pi=Pi+ piT, (17) 
for i = 1,2. Note that, unless specified otherwise, the variables are functions of t. Let 

6 = \ql (18«) 
£2 = ^ 2 2 , (186) 
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[18d) 



Here 

H= -6 -6+6+6 (19) 
is a linear function of the £s. Differentiating the £s with respect to time, we get 

6 = 9i9i, (20 a) 

6 = 9292, (206) 

6 = ©92 = 92 (P2 -9i), (20c) 
6 = 6+6- 6, (20d) 
on making use of (0), together with the conservation of H. The corrector is given by 

(21) 



T 



6(* + r)=6 + -(6 + 6), 

for z = 1, ... ,4, where 6 is simply ( |18q| ) evaluated at p« and £ + r. Inverting, the 
new values of % and pj can be expressed in terms of 6 as 



91 = sgn(9i)v/2^, 

92 = sgn(g 2 )v / 26> 



and, on using (|15q ) and (|156|) , 



Pi 



2(1 -/x) 2/i 
-9a + sgn(pi + 92)4/26 + — + — , 



r-2 



P2 = 9i + sgn(p 2 - 9i) V26- 



(22a) 
(226) 

(23 a) 
(236) 



This example assumes that the mass of one body is negligible to the other two 
masses and that the other two masses are travelling in circular orbits. The rest of this 
paper discusses the general case of three or more bodies: no restrictions are placed on 
the masses of the bodies, and their orbits do not have to be circular, or even periodic. 



5. General Three-Body Problem 



The derivation of the equations of motion of the general three-body problem in a plane 
is described in Refs. [|13|], |T3]1 , and |i~5 |. 

Given three bodies m 1; m 2 , and with position vectors r 1; r 2 , and r 3 , where 



each Vi is at location 



V 



Xi,yi), define r# = 
Gm\m2 Gm 2 m3 



rj - n, for i,j 
Gm 1 777.3 



where G is the gravitational constant and r»j 
between the ith and jth bodies. 



ri3 



1,2,3. The potential is 

(24) 

a/ (xj — Xi) 2 + (yj — yi) 2 is the distance 
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The system consists of three second-order differential equations, 
dV Gm\m 2 (r 2 — T\) Gm\m 3 (r 3 — 7*1) 



niiri = — 



m 2 r 2 



dri 
dV 



12 



13 



Gm x m 2 (r x - r 2 ) Gm 2 m 3 (r 3 - r 2 ) 



m 3 r 3 



dr 2 

dV 
cV 3 



+ 



21 



23 



Gmim^ri - 7*3) Gm 2 m 3 (r 2 - r 3 



+ 



(25 a) 
(25b) 
(25 c) 



32 



3 '31 

These equations conserve the total linear momentum Y2i=i m %^% (which allows us to 
fix the center of mass at the origin) and total angular momentum X^=i r i^ m &i- The 
Hamiltonian 

3 



(26) 



i=i 



where V is given by (p4[), is also conserved. We exploit the constancy of the linear 
momentum and center of mass position to reduce the number of degrees of freedom 
in the problem. It is convenient to implement this reduction by converting to Jacobi 
coordinates (e.g., see Refs. [16|], |l7j, and ||18|| ). The remaining constraints of constant 
total angular momentum and energy are built into the conservative integrator by 
transforming to a frame where these invariants are linear. 

Letting r = r 2 — r\ = (r x ,r y ), M = mi + m 2 + m 3 , and fi = m\ + m 2 , the 
location of the center of mass of mi and m 2 is seen to be at /i _1 (miri + m 2 r 2 ), or, since 
mir 1 + m 2 r 2 + m 3 r 3 = 0, at — /i _1 m 3 r 3 . Let p = (p x , p y ) be the vector from the center 



H 1 m z r 3 



of mass of the first two bodies to the third body. Then p = r 3 
and we find 

r 2 - n = r, 

r 3 - 7*1 = p + m 2 p~ 1 r, 

r3 - r 2 = p - m\pT x r. 

In these coordinates, following (^6|) , the Hamiltonian can be written as 

H = \gx(rl + rj) + -g 2 (pl +p 2 y ) + V 

in terms of the reduced masses g\ = m 1 m 2| u" 1 and g 2 = m 3 p/M, where V is given 
by ©■ 

Define r x = rcos6, r y = rsin^, p x = pcosG, and p y = psinB. In these polar 
coordinates, the Hamiltonian can be rewritten 

.2 p2 £ 2 L 2 



Mp- X r 3 

(27a) 
(27b) 
(27c) 

(28) 



+ 



+ v(r, P ,e,e), 



(29) 



2 9l 2g 2 2 9l ri 2g 2 p* 
where p is the linear momentum of the first reduced mass, I is the angular momentum of 
the first reduced mass, P is the linear momentum of the second reduced mass, L is the 
angular momentum of the second reduced mass, and V = V(r, p,9,Q) is the potential 
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energy of the system. The Hamiltonian H and the total angular momentum i + L are 
conserved, and the center of mass remains at the origin for all time. 
The equations of motion in polar coordinates are 
. OH p . 8H I 

8p g\ at g\r l 

. _8H e _8V I- _ dH _ dV ( 30 M 
^ dr g±r 3 dr ' 89 86 ' 

8H P ■ 8H L 

8P g 2 8L g 2 p 2 

■ _ 8H L? 8V ■ 8H 8V 
5.1. Integration 

The variables can be transformed as 

p 2 £ 2 P 2 L 2 

2g l 2g x r 2 2g 2 2g 2 p 2 

& = v, ^ = P , & = & = l, t 7 = e, 6 = e,(3i6) 

so that the conserved Hamiltonian becomes a linear function of the transformed 
variables: H = ^ + £ 2 + £3- The time derivatives become 

6 = 7 + (32a) 

& = — + fc^, (32.) 
92 g 2 p A 

■ 8V 8V ■ 8V 8V ■ , , 

6 =7F r + 7» s + 7*r' , + 5e e ' < 32c) 
6, = p, & = & = i, & = & = e. (32<i) 

The integration procedure is an extension of the method used for the Kepler problem. 
We can invert to find the original variables as follows, 

p = U £ = £ 5 , L = £ 6 , 9 = &, = £ 8 , (33a) 
r = g(^,p,e,Q), (336) 



p = sgn(p)^/2 5 i(ei-2^2 ). ( 33c ) 



P = sgn(P)^/2^ 2 (6-^). (33c?) 



An Exactly Conservative Integrator for the n-Body Problem 

1 



10 



* o 




Figure 1. Conservative predictor-corrector solution for a general three-body 
choreography. 



The value of the inverse function g defined by V(g(^, P, 0, ©), p, 0, ©) = £3 is determined 
at fixed p, 8, O by Newton-Raphson iteration, using the predicted value f as an initial 
guess. 

In Fig. [I] we use our conservative predictor-corrector to illustrate the remarkable 
three-body figure-eight choreography discovered by Chenciner and R. Montgomery (|TJJ 
and located numerically by Simo ||20|| . The dots indicate the initial positions of the three 
unit masses. The gravitational constant G is taken to be unity and the initial conditions 
are those cited in |1| : 

ri = (0.97000436, -0.24308753), r 2 = (0, 0), 
r 3 = (-0.97000436,0.24308753). 

fi = (0.46620369, 0.43236573), r 2 = (-0.93240737, -0.86473146), 

r 3 = (0.46620369,0.43236573). (34) 

We used a fixed time step of r = 10 -4 and integrated for a complete choreographic 
period, during which each mass travels once around the figure eight. 



6. General n-Body Problem 



The Jacobi coordinates can be extended to n > 2 bodies in a plane [18, 16]. Let each 



body of mass m; have radius vectors r,, where i — 1, . . . , n. Define r%j = rj — Vi as the 
vector joining rrii to rrij. Also define Ci to be the center of mass of the first i bodies, 
where i — 2, . . . , n, and choose the origin of the coordinate system so that C n = 0. Let 
the vectors pi be defined such that 



P2 = r*i2, 



(35 a) 
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p 3 = r 3 -C 2 , (35b) 

(35c) 

P» = »n-C n _i. (35d) 



Also, 



£-1 



where 1 < k < £ < n, and Mj = 5^ =1 mfc.§] 

The reduced masses are 
m 2 m 1 

92 = -Ti—, (37) 



M 2 



m 3 (m 2 + mi) 
#3 = » 



(39) 

= (40) 

The equations of motion in polar coordinates are just an extension of the three-body 
problem: 

a£ Vr (41a) 

* = W- = A> (41^) 
Pi = ~— = :r3 - — > ( 41c ) 

(41*0 

where pj, 0j, and ^ are the radius, angle, linear momentum, and angular momentum, 
respectively, of the ith reduced mass, for i — 2, . . . , n. The potential is defined to be 

V = -J2 (42) 

i,i=i 

and the total kinetic energy is 





9iPi 




dH 




dV 


dpi 


9iPf 


dpi' 


dH 


dV 




09, 







It is easy to verify that the Hamiltonian H = K + V is conserved by ( [41a ). The total 



angular momentum X^=2 ^ ^ s a ^ so conserved, and the center of mass remains at the 
origin for all time. 

§ Here p\ is a dummy variable that cancels out in the expression for m- 
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Transform (p,6,p,i) to (£,0, 77,-6), where 

C2 = V, (44a) 
Q = pi, for % = 3, . . . , n, (446) 

*L + J? 
2& 2<^ 2 

Note that if is a linear function of the transformed variables: 



Vi = + for i = 2,...,n. (44c) 



# = 5> + C2, (45) 

as is the total angular momentum L = J2i =2 £i- The time derivatives of C, and 77 are 
given by 

• >A fdV . dV ■ \ 

Ci = pi, for i = 3, . . .,n, (466) 
ft = gft + <!gfrz»g* for i = 2,...,„. (46c) 

& #iPi 

The predictor equations are 

Pi = Pi + PiT, #i = ®i + ^t, (47a) 
Pi = Pi + PiT, £i = £i + £iT (476) 

and the corrector is given by 

(i(t + t) = Ci + ^(6 + + r) = 0i + ^ + h), (48a) 

7^ + t) =^ + 1(^ + 77-), £i{t + T) =ei + ^(£i+'ii), (486) 
for i = 2, . . . , n. 

One then inverts to get the original variables as functions of the temporary 
transformed variables: 

Pi = Ci for i = 3, ... ,n, (49a) 
P2 = g(C2,P3,---,Pn,0), (496) 



<?2 



^ = sgn(p;)^/2^ ^r/i - ^2 J , for i = 2, . . . , n. (49c) 

The value of the inverse function g defined by 

V(g(C2, P3, • • • , Pn, 0),p 3 ,..., p n , 6) = ( 2 (50) 
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is determined at fixed p 3 , . . . , p n , with a Newton-Raphson method, using the predicted 
value p2 as an initial guess. 

In Fig. 0, we illustrate Simo's four-body choreography [20|. The motions of one 
of the four unit masses as determined by the predictor-corrector and conservative 
predictor-corrector algorithms are compared, using the fixed time step r = to 
integrate the system from time t — to t — 11.5. The gravitational constant G is taken 
to be unity and the initial conditions are given by 

ri = (1.382857, 0), r 2 = (0, 0.157030), 
r 3 = (-1.382857, 0), r 4 = (0, -0.157030), 
ri = (0, 0.584873), r 2 = (1.871935, 0), 

r 3 = (0, -0.584873), r 4 = (-1.871935, 0). (51) 

As r is decreased, the predictor-corrector solution converges to the conservative 
predictor-corrector solution obtained with a large time step. This emphasizes that the 
conservative predictor-corrector can be viewed as a finite-time-step generalization of the 
conventional predictor-corrector, as argued in Ref . |IJ . We also compare these solutions 
to a symplectic map based on the simple second-order kinetic-potential energy splitting 

Pi=Pi~ \^ V ^ ?2, • • • , Qn), 

d 

q'i = qi + T—K(p 1 ,p 2 , . . .,p N ), 

Pi -m- (52) 

to evolve the canonical variables (qi,Pi) to (q'^Pi), for i = 1,...,N. This second- 
order scheme, which is implemented as the method SKP using Varadi's NBI code with 
t = 10 -3 , is similar to the one described by Forest and Ruth |7|, §], with the roles of the 
coordinates and momenta interchanged. 

In Figure |3], we compare the root-mean-square error in the computed trajectory 
between t = and t = An (twice the period of the choreography), for each of these 
integration algorithms. The error was computed relative to a fifth-order Runge-Kutta 
integrator with time step r = 10~ 5 . Of the three other solutions, we note that the 
conservative predictor-corrector trajectory is the most accurate. For general n-body 
integrations, our conservative algorithm was also observed to be more accurate than the 
second-order symplectic Wisdom-Holman scheme ^TQ [22], 0] , but this is expected since 
the latter applies only to small perturbations of Keplerian orbits. 



7. Conclusion 



Conservative integration algorithms can reduce the computational effort required to 
integrate a system of ordinary differential equations. For example, when the total 
energy and angular momentum of the n-body problem is conserved, it is possible to 
obtain accurate trajectories using a larger time step than with conventional integration 
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Figure 2. The predictor-corrector (dotted line), symplectic (dashed line) and 
conservative predictor-corrector (solid line) solutions for a four-body choreography. 




Figure 3. Root-mean-square error of the predictor-corrector (dotted line), symplectic 
(dashed line) and conservative predictor-corrector (solid line) solutions in Fig. ^. 



methods. This is particularly relevant for extremely long-time integrations. In contrast, 
symplectic methods typically predict a total energy that oscillates about the correct 
value. In some cases, these oscillations can eventually lead to large excursions from the 
mean value, similar to random walk diffusion. In addition, there are certain statistical 
mechanical systems (such as equipartition states of inviscid fluids) where the final mean 
state is a function of only the initial values of the invariants; for these systems, a 
conservative integrator is clearly preferable to a symplectic algorithm. However, as 
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these methods put the integration error into different places, the integration method 
that is most suitable for a given system ultimately depends on the nature of the 
physical problem, the integration time scale, and the kinds of questions addressed by 
the numerical simulation. 

In the case of the n-body problem for planar motion, there are six invariants, all 
of which need to be considered during the integration. Jacobi coordinates were used 
to reduce the system to an (n — l)-body problem in which the linear momentum and 
center of mass constraints are implicitly built in, leaving fewer conservation laws to be 
explicitly built into the algorithm. In Jacobi coordinates, the kinetic energy term of the 
Hamiltonian remains in diagonal form (a sum of squares); this makes it easy to express 
the Hamiltonian as a linear function of new variables. 

Future work in this area should include extending the numerical code to the full 
three-dimensional case and regularizing the potential terms to handle collisions and 
close approaches. One could also build in precession, nutation, and tidal effects into the 
equations of motion. 

This work was supported by the Natural Sciences and Engineering Research Council 
of Canada. 



Appendix A. Error Analysis 



Here we describe the local error analysis of the second-order conservative predictor- 
corrector scheme given by (|2~q ) and ([|). We assume that both / and T are analytic 
functions and that the points where T' vanishes are isolated. For notational simplicity, 
we restrict the analysis to the autonomous one-dimensional ordinary differential equation 
dx/dt = f(x), for which the exact solution is given by 



2 3 

x (t + t)=x + r/fao) + ^f{x )f{x ) + T-f( Xo )f( XQ ) + C(r 4 ) 

2 

The conservative predictor-corrector scheme 
x = x + rf(x ), 

£(t + r) = £0 + T - [T'(x )f(x ) + T'(x)f(x)\ 



yields the solution 



T 



2 f2 



T f f + T'f + (T f f)Tf + (T'f) 



* r'f 



0(t 3 



£0 + rT'f + — {T"f 2 + T'f'f) 



T 

T 



(T"'f + 2T"/'/ 2 + T'f'f) + 0(r 4 



(A.l) 

(A.2) 
(A.3) 



(A.4) 



where the expressions on the right-hand side are all evaluated at xq. The new value of 
x is then given by 

x(t + T)=T- 1 ^(t + r)) 
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.3 



rT'f+ T -{T"f 2 + T'f'f) 



— (T"'f + 2T"ff + T'ff) + 0{t a 



rT'f + ^-(T"f 2 + T'f'f) + 0(r 3 



+ V lw (£ ) [rT7 + 0(r 2 )] 3 + 0(r 4 



By implicitly differentiating the identity T 1 (T(a;)) = s, it follows that 

-| /TV/ nrpll2 rplll 

rp-llllp \ _ __ rp-llllfp \ _ 

rp,1 ± ISO; rp /3 , J- \^0J rp /5 rp^, 



rp — lf 



(ft 



(A.5) 



(A.6) 



so that ( |A.5| ) simplifies to 



x{t + r) = x Q + r/(x ) + -z-f'(x Q )f(x ) 



+ 



/"(xo)/ 2 (x ) + 



T'"(x ) 



/ 3 M 



(A.7) 



3T'(x )" 

By setting T(x) = x, we obtain the usual error estimate for the conventional predictor- 
corrector. We see that both conservative and conventional schemes are accurate to 
second order in r; moreover, for quadratic transformations like T(x) = x 2 , which in 
light of Lemma [I] are often useful for enforcing energy conservation, the conservative 
and conventional schemes agree through third order in r. The appearance of T'(xq) in 
the denominator of the third-order (error) term emphasizes that we must exercise care 
at singular points of T. Near these points, either a conventional scheme can be used 
or the time step can be reduced, as previously remarked 0. Should T'(xq) = 0, the 
estimate ( |A.7| ) should be replaced by 



x(t + 



T 



-i 



T(x ) + -T'(x)f(x] 



(Ai 



which is guaranteed to have a solution for sufficiently small r if the points at which T 1 
vanishes are isolated. The transformation is then invertible at x(t + r), allowing the 
integration to be continued beyond the point of singularity. 
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